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ABSTRACT:  FNOL3  Is  a  Fortran  IV  subprogram  which  uses  fourth  order  Runge-Kutta 
and  Adams-Moulton  methods  to  solve  up  to  30  ordinary  differential  equations  with 
initial  conditions.  Options  allow  control  of  error,  step-size,  print  frequency  and 
integration  method.  There  is  a  discussion  on  the  integration  methods  and  their 
error  terms.  Listings  of  the  subprogram  and  sample  problems  and  their  results 
as  run  on  NOL's  CDC  6400  are  included. 
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FNOL3,  A  COMPUTER  PROGRAM 
TO  SOLVE  ORDINARY  DIFFERENTIAL  EQUATIONS 


I.  INTRODUCTION 

This  report  describes  a  computer  program  FNOL3  for  the  numerical  integration 
of  ordinary  differential  equations  with  initial  conditions.  These  equations  are 
reduced  by  the  user  to  a  system  of  first  order,  simultaneous,  ordinary  differential 
equations  with  initial  conditions  and  are  solved  by  using  fourth  order  Runge-Kutta 
and  fourth  order  Adams-Moulton  predictor-corrector  methods.  FN0L3  is  the  successor 
to  the  Naval  Ordnance  Laboratory  ordinary  differential  equation  solver  FN0L2 
and  uses  many  features  familiar  to  users  of  FN0L2  (the  differences  are  provided 
in  Appendix  E) .  See  reference  (a) . 

FN0L3  is  formulated  so  as  to  be  quite  separate  from  any  particular  application. 
Options  make  available  a  convenient  flexible  package  that  can  be  used  whenever  the 
problem  is  expressible  as  a  system  of  equations  as  described  above. 

The  program  is  written  in  the  FORTRAN  IV  language  for  the  operating  system 
currently  used  on  the  Laboratory* s  CDC  6400  computer. 

II.  DESCRIPTION  OF  THE  METHODS 

The  discussions  that  follow  in  this  section  and  in  Section  III  are  taken 
from  reference  (b) .  The  notation  used  omits  commas  from  subscripts  (f^  )  except 
when  there  might  be  an  ambiguity  (f  ^  n+1)  •  Y  is  used  for  true  values  of  the 


dependent  variable  and  y  for  computed  values. 


Let  the  system  of  equations  to  be  solved  be  given  in  the  form: 


_  dy± 


"  fi(x,y1,y2,...,yN) 


i  =  1,2,. ...N 


(la) 


yi(x0)  *  yi0 


1 


N 


(lb) 


T  -  T(x,y1,y2,...,yN,y'1,y,2, 


•  •  •  , 


(lc) 
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where  T  is  the  termination  condition  and  y^o’^201 


yN0 


are  the  initial 


conditions  at  x  =  x^. 

Let  y.  be  the  value  of  y.  and  f .  be  the  value  of  f .  at  x=x  .  Let  h  be 
in  J  i  in  in 

the  step  size  in  the  independent  variable  x.  The  Runge-Kutta  method  uses  the 
following  formulae  with  the  appropriate  initial  conditions  to  go  from  step  n  to 
step  n+1. 

kU  ■  £i(Vyln . yN„>  1 

ki2  ’  . 1 

ki3  ‘  . 1 

ki4  ’  fi(Vh>yln+hk13 . yN„+hkN3)  1 

yl,o+l  "  yin+6(kil+2k12+2k13+k14)  1 

The  Adams-Moulton  predictor-corrector  method  uses  the  following  formulae 

to  compute  the  values  of  y,  using  the  values  of  y.  «,  y.  _ ,  y.  ,  ,  and 

l,n*rl  i,n-o  i  ,n-z  i,n-l 

yin’ 

y,(p)  =  y,  +^r(55f,  -59f .  ,+37f , 

Ji,n+1  J±n  24  in  i,n-l  i. 


:  x=x  .  Let 
n 

h  be 

method  uses 

the 

go  from  step 

n  to 

1,2,.  .  .  ,N 

(2a) 

1,2 . N 

(2b) 

1,2,.  . .  ,N 

(2c) 

1,2,...  ,N 

(2d) 

55 

CM 

r— 1 

(2e) 

i,n+l 


y  +-^(9f(P) 

yin  24v*  i,n 


,.+19f,  -5f, 
n+1  in  i 


,+f 


-9fi,n-3> 

i  =  1,2,. ,.,N 

(3a) 

i.n-25 

i  =  1,2,. . . ,N 

(3b) 

where  the  values  of  y^  at  x^,  x^ ,  and  x^  are  found  by  using  the  Runge-Kutta 

formulae,  is  the  first  point  at  which  the  Adams-Moulton  method  is  used. 

Observe  that  the  Adams-Moulton  formulae  only  require  two  derivative 
evaluations  to  go  to  the  next  step.  The  Runge-Kutta  method  requires  four  deri¬ 
vative  evaluations  for  each  step. 


Ill .  THE  ERROR  TERMS 

We  will  assume  in  this  section  only  one  ordinary  differential  equation,  so 

the  i  from  y.  and  f .  is  omitted.  The  error  term  for  the  Adams-Moulton  method 
Jin  in 

is  as  follows : 
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y(^  ■  y  +fe<55f  ~59f  i+37f  o~9f  J 

n+l  24  n  n-1  n-2  n-3 

y(ii  =  y  +57<9f(^+19f  ~5f  !+f  o) 

n+l  7n  24'  n+1  n  n-1  n-2 


(la) 

(lb) 


:  (p) 


where  p  stands  for  predicted  value  and  c  for  corrected  value,  f  ‘  ,  means 

n+l 

(P) 


calculate  f  at  using  y^£' 

If  the  calculated  values’  of  yi,y2,'’’’yn  anc*  accordingly,  of  f  ]_  >f  2  » •  •  •  *^n 
were  exactly  correct,  then  the  true  ordinate  at  xn+i>  say  Y^+^  would  satisfy  the 


equations 


Vl  ■  V2i<55£„-59£n-l+37f„-2-9£„-3)+27i  1,5 


(2a) 


Vl  ■  >’n+2Ji<9£„+l+19£„-5£n-l+£n-2)  ‘  #  ^  *  V 


(2b) 


where  and  both  lie  between  x  ^  and  x  y  (£-)  is  the  fifth  derivative 

i  l  n-J  n+l  1 

of  y  at  x 


So  it  follows  from  (1)  and  (2)  that 

r(p) 


v  _  251  ,5  V,c  N 

Yn+1  yn+l  720  h  y  (C1) 


(3a) 

(3b) 


Y  -y^  *  — (f  -f  ^)  — —  yV(£  ) 

n+l  yn+l  8^  n+l  n+l'  720  y  U2; 

Let  F(xn+1*Yn+1)  =  fn+l^Xn+l,Yn+l^ ’  aPP1yin8  the  law  of  the  mean 

fn+l  ~  fn+l  =  F(xn+l,Yn+P  "  F(xn+l,Yn+l)  =  (Yn+l‘yn+l)Fy(Xn+l’rin+l)  (4) 

where  ri  , ,  is  between  Y  ,,  and  y^?,. 
n+l  n+l  n+l 

It  is  assumed  that  h  is  sufficiently  small  to  ensure  that 


8  lFy*Xn+l’nn+l^  <<  1 


(5) 


and  also  that  yV(x)  does  not  vary  strongly  for  x^_^  <  x  <  xn+j  so  that  y^(£^) 
and  can  be  equated.  Equations  (3)  lead  to  the  useful  approximate  relation 

I20(Y  -y(p))  =  h5  yV(5) 

251  n+l  yn+l;  7 


n+l 


(c)  =  3h/y  (P)%  f  /x  _  (Z2£wy  -v(p\ 

8(n+l  Yn+l'  y(  n+l’nn+l7  (720H251Mn+l  yn+l' 


n+l 


Yn+1  yh+l  _(38  Fy(xn+l’nn+l)  251J (Yn+l  yn+^ 
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so  from  (5) 


n+1 


Yn+1  + 


.  v(0  ~  19 


n+1 


(Y  -y(p) ) 
251^  n+1  yn+l' 


-Aly  y(0  +  JL9  (p) 

251  n+1  ~  yn+l  251  yn+l 


270  ~  270  (c)  19  (c)  ,  19  (p) 

251  n+1  251  yn+l  “  251  yn+l  251  yn+l 


Therefore, 


Y  -  y^  +  —  (y^  -  y^) 
n+1  ~  yn+l  14  vyn+l  yn+l' 


(6) 


This  is  the  equation  programmed  when  the  Adams -Moulton  method  is  selected. 
The  fourth  order  Runge-Kutta  formulae  are 


f(x  ,y  ) 

(7a) 

n  n 

f (x  4%h,y  +^hk. ) 
n  'n  1 

(7b) 

f(x  +^h,y  44shk  ) 
n  n  z 

(7c) 

f(x  +h,y  +hk») 
n  n  3 

(7d) 

‘  yn  +  5«‘l+2kJ+2k3*V  ' 

(7e) 

No  simple  expressions  are  known  for  the  precise  trincation  errors  in  the 

preceding  formulae.  An  estimate  of  the  error  can  be  obtained,  in  practice, 

in  the  following  way.  Let  the  truncation  error  associated  with  a  formula  of 

4th-order  accuracy  in  advancing  from  x^  to  that  at  xr+^  =  x^  +  h  in  a  single 

step  be  denoted  by  c^  h"\  Also  suppose  that  c^  varies  "slowly"  with  n  and  is 

nearly  independent  of  h  when  h  is  small.  Then  if  the  true  value  of  y  at  xn+^ 

is  denoted  by  Y  j.,  the  value  obtained  by  two  steps  starting  at  xn_^  by 

(2h) 

and  the  value  obtained  by  a  single  step  with  doubled  spacing  2h  by  there 

follows  approximately 
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Y  -  y0^  S  2c  h5 

n+1  •'n+l  n 

Y  -  y(J?)  =  c  (2h) 5 

n+1  ■'n+l  n 


(8a) 

(8b) 


when  h  is  small. 


From  equation  (8b) 

Y  -  y(^°  S  25c  h5 
n+1  Jn+1  n 


so 


i  -  y^2h^  =  2^(Y  -y^) 

n+1  yn+l  Un+1  yn+lJ 


Therefore , 


Vl  S  y;+l  + 


,  %  /w.y(2h> 

(h)  yn+l  yn+l 


15 


(9) 


This  is  the  equation  programmed  when  the  Runge-Kutta  method  is  selected. 


IV.  AUTOMATIC  ADJUSTMENT  OF  STEP  SIZE 

FN0L3  has  the  option  of  automatically  varying  the  step  size  h,  to  hold 
the  truncation  errors  within  bounds  fixed  by  the  user.  The  absolute  truncation 


errors  for  the  Adams-Moulton  method  from  III-6  are 
y(p)  _v(c) 

AE. 


yi,n+l  yi,n+l 


14 

and  the  relative  truncation  errors  are 

AE . 


RE, 


,<c) 

i,n+l 


1,2, . . .  ,N 


i  =  1,2, ...  ,N 


The  absolute  truncation  errors  for  the  Runge-Kutta  method  from  III-9  are 

i  =  1,2,. ...N 


>>  _y<2h> 

yi,n+l  yi,n+l 

“i  "  15 


and  the  relative  truncation  errors  are 

AE. 


RE±  = 


y(h) 

yi,n+l 


i  =  1,2, .  ..  ,N 


(1) 


(2) 


(3) 


(4) 
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To  determine  the  step  size  at  each  step,  let  E^.Ej , . . . ,E^  be  the  errors 

whether  they  are  absolute  or  relative  errors  and  let  XNE  be  a  step  size  control 

— YNF— ^  i  i  _Y"NF 

parameter  (see  page  8).  Then  if  10  <  E^|  <  10  for  i  =  1,2,...,N, 

i  i  -XNE-3 

the  step  size  is  not  changed.  But  if  |E^|  <10  for  all  i,  then  let 


HB  -  minimum 


10 


-XNE-1.5 


>  1 


i  IeJ+10-11 

and  the  step  size  is  increased  to  HB^r-  h.  If  |Ej|  >  10-™  for  some  j  then 


(5) 


HB  =  minimum 


10 


-XNE-1.5 


<  1 


J  lEjl+10 


and  the  step  size  is  decreased  to 


-11 

HB^r<  h. 


(6) 
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V.  PROGRAMMING 

The  user  must  write  a  calling  program  hereafter  called  MAIN,  and  three 
auxiliary  subprograms.  The  latter  are  usually  called  DERIV,  TERM,  and  OUTPUT 
and  are  described  later  in  this  section.  Besides  calling  on  FN0L3  and  providing 
it  with  initial  values,  MAIN  must  contain  an  EXTERNAL  card  which  denotes  that 
DERIV,  TERM  and  OUTPUT  are  the  names  of  subprograms,  not  variables. 

The  calling  sequence  for  FN0L3  is 

EXTERNAL  DERIV,  TERM,  OUTPUT 

CALL  FNOL3  (J ,N ,G,L ,M,XNE ,X,Y ,D , DERIV .TERM, OUTPUT) 

J:  (INPUT, INTEGER) 

This  parameter  indicates  the  integration  method. 

J  =  1  Use  Runge-Kutta  method  of  integration  to  termination.  Truncation 
errors  are  not  calculated;  the  step  size  G  is  not  adjustable. 

J  =  2  Use  Runge-Kutta  for  the  first  three  steps,  then  Adams-Moulton 
for  the  remainder  of  the  interval  of  integration.  Truncation  errors 
are  calculated.  The  step  size  is  adjustable  unless  XNE  =  0.  If  the 
step  size  is  adjusted,  new  starting  values  are  obtained  through  the 
Runge-Kutta  method. 

J  =  3  Use  Runge-Kutta  throughout.  The  truncation  errors  are  calculated; 
the  step  size  is  adjustable  unless  XNE  =  0. 

N:  (INPUT, INTEGER) 

This  is  the  number  of  simultaneous  first  order  differential  equations 
to  be  solved.  The  maximum  number  of  equations  is  30. 

G:  (INPUT, REAL) 

This  is  the  initial  step  size;  upon  return  from  FN0L3,  G  retains  its 
original  value. 
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L:  (INPUT, INTEGER) 

This  is  the  number  of  Y's  beyond  (N+3)  to  be  written  in  the  routine 
OUTPUT.  These  additional  Y's  should  be  calculated  in  the  routine 
TERM,  beginning  with  Y(N+3+l)  to  Y(N+3+L) . 

M:  (INPUT, INTEGER) 

This  is  the  number  of  accepted  steps  taken  between  calls  to  routine 

OUTPUT.  If  M  =  0,  then  printing  is  determined  by  values  assigned  to 

Y(N+1)  and  Y (N+2) .  Y(N+1)  =  f(X,Y(l) . Y (N) ,D(1) , . . . ,D (N) )  and  is 

defined  in  routine  TERM.  To  generate  a  call  to  OUTPUT,  Y(N+1)  must 
change  by  an  amount  greater  than  Y(N+2)  since  the  last  call  to  OUTPUT. 
Y(N+2)  is  assigned  a  constant  value  in  the  routine  which  calls  FN0L3. 

XNE:  (INPUT, REAL) 


This  is  the  step  size  control.  The  step  size  is  unchanged  if  the  worst 

-XNE -3  —XNE 

of  all  the  errors  lies  within  the  window  [10  ,10  ].  The  step 

—XNE— 3 

size  is  increased  if  the  errors  are  all  less  than  10  .  The  step 


size  is  decreased  if  for  some  differential  equation  the  error  is  greater 

in-XNE 
than  10 


If  Y(N+3)  <  0.  and  XNE  4  0.,  the  automatic  adjustment  of  the  step 
size  is  a  function  of  the  absolute  errors. 

If  Y(N+3)  =  0.  and  XNE  4  0.,  the  automatic  adjustment  of  the  step 
size  is  a  function  of  the  relative  errors. 

If  Y(N+3)  =  e  >  0.  and  XNE  4  0.,  the  automatic  adjustment  of  the 
step  size  is  a  function  of  the  relative  errors  where  the  relative 
errors  are  equal  to  the  absolute  errors  divided  by  the  maximum 
(Y(N+3),  |Y(I) | ) .  This  option  removes  the  possibility  of  using 
"small"  functional  values  to  compute  relative  error,  otherwise  this  option 
is  identical  to  the  previous  option  and  is  to  be  preferred  over  it. 
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If  XNE  =  0.,  the  step  size  G  is  not  adjustable. 

X:  (INPUT, OUTPUT, REAL) 

This  is  the  independent  variable.  An  initial  value  must  be  specified 
before  calling  FN0L3.  FN0L3  will  return  with  X,the  terminal  value  of 
the  independent  variable. 

Y:  (INPUT, OUTPUT, REAL) 

This  is  the  name  given  the  solution  array.  Y  must  be  dimensioned  at 
least  Y(N+3+L).  Initial  values  for  Y(l) ,Y(2) , . . . ,Y(N)  must  be 
specified  before  calling  FN0L3.  If  L  >  0,  then  Y(N+3+l) , . . . ,Y(N+3+L) 
may  be  used  to  calculate  additional  values  in  the  routine  TERM.  Upon 
returning  from  FN0L3 ,  Y(l) , . . . ,Y(N) ,Y(N+3+l) , . . . ,Y(N+3fL)  have  the 
values  computed  at  the  terminal  value  of  X. 

D:  (OUTPUT, REAL) 

This  is  the  name  given  to  the  array  where  the  derivatives  are  stored 
and  must  be  defined  in  routine  DERIV.  D  should  be  dimensioned  D(N) . 
Upon  returning  from  FN0L3,  D(l) , . . . ,D(N)  have  the  values  computed  at 
the  terminal  value  of  X. 

DERIV: 

In  this  routine  the  user  must  compute  the  N  derivatives.  The  general 
form  is 

SUBROUTINE  DERIV  (X,Y,D) 

DIMENSION  Y (1) ,D(1) 

D(l)  =  ••• 

D(2)  =  • •  • 

D(N)  =  • • • 

RETURN 

END 
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No  D  beyond  D(N)  should  be  calculated.  If  desired,  additional  data  may 
be  passed  between  this  routine  and  the  other  user  written  routines  via 
COMMON  statements. 

TERM: 

The  user  evaluates  the  termination  criterion  (T  =  T(X,Y(1) , . . . ,Y(N) , 

D(l) , . . . ,D(N) )  in  this  routine.  Auxiliary  values  such  as  Y(N+1)  and 

those  required  for  plotting  purposes  should  be  calculated  here  because 

this  routine  is  entered  only  once  per  accepted  step  before  termination. 

A  termination  loop  of  at  most  four  iterations  starts  when  T  undergoes 

a  change  of  sign.  After  the  very  first  step  of  the  integration, 

termination  may  occur  without  looping  when  | T !  <10 

The  general  form  is 

SUBROUTINE  TERM  (X,Y,D,T) 

DIMENSION  Y(l) ,D(1) 

•J*  sr  •  •  • 

RETURN 

END 

OUTPUT: 

This  user  routine  is  entered  at  the  beginning  and  end  of  the  complete 

integration.  If  M  ^  0,  it  is  entered  every  Mth  accepted  step.  If 

M  =  0  then  Y(N+1)  and  Y(N+2),  discussed  under  parameter  M  determine 

the  print  frequency. 

The  general  form  is 

SUBROUTINE  OUTPUT  ( X,Y,D, ERROR, N,L,H) 

DIMENSION  Y(l) ,D(1) ,ERR0R(1) 

PRINT 

RETURN 

END 

X,Y,D,N,  and  L  are  the  same  as  in  the  calling  sequence  of  FN0L3. 

ERROR  is  the  name  given  the  array  which  contains  the  absolute  errors. 
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H  is  the  step  size  used  to  get  to  the  present  step.  If  the  magnitude  of  a 
new  step  size  becomes  less  than  «  10  ^  *  |x|  +  10  then  two 
steps  are  taken  using  the  Runge-Kutta  method  with  step  size  equal  to 
or  -H^  depending  on  the  direction  of  integration.  If  the  errors 
are  too  large,  the  two  steps  will  be  accepted  and  two  more  steps  are 
taken  using  the  Runge-Kutta  method  with  a  step  size  computed  in  the 
same  way  as  above  using  the  current  x.  If  the  errors  are  still  too 
large  after  the  above  procedure  is  done  ten  times,  then  an  error 
message  is  written  and  the  program  stops  after  going  to  OUTPUT.  If 
at  any  step  the  errors  are  not  too  large,  then  FN0L3  continues  in  the 
normal  manner  after  resetting  the  counter  back  to  zero. 
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VI.  EXAMPLES 

The  examples  section  of  this  report  contains  three  problems  which  illustrate 
the  available  options  of  FN0L3. 


Example  No.  1 

Solve  the  first  order  linear  differential  equation 

=  uA  cos(ux) 

in  the  interval  [0,3tt]  for  the  initial  condition 
y(0)  =  0 

Choose  the  case  u  =  A  =  1. 

Use  the  three  available  methods  of  FN0L3  to  numerically  integrate  the  given 
differential  equation  over  [0,3ir]. 

(I)  Use  the  fourth  order  Runge-Kutta  method  with  a  constant  step  size 

of  .1. 

(II)  Use  the  fourth  order  Adams-Moulton  predictor-corrector  formulae 
with  an  adjustable  step  size  to  hold  the  relative  truncation  error  in  the  window 
[10"8,10"5]. 

(III)  Use  the  fourth  order  Runge-Kutta  method  with  an  adjustable  step 

—8  —5 

size  to  hold  the  relative  truncation  error  in  the  window  [10  ,10  ]. 

The  analytic  solution  of  the  given  equation  with  respect  to  the  given  initial 
condition  is 

y(x)  =  A  sin(u>x). 

The  actual  value  of  sin(x)  is  computed  for  each  x  and  compared  to  the  calculated 
value  of  the  solution.  At  every  fifth  accepted  step,  the  solution,  sin(x) ,  and 
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the  difference  of  the  actual  and  calculated  solutions  are  printed.  These  results 
are  in  Table  1. 

The  main  (calling)  program  and  the  necessary  subroutines  to  accomplish  the 
solution  of  the  stated  problem  in  the  manner  indicated  are  in  Appendix  A. 

Example  No.  2 

Solve  the  second  order  linear  differential  equation 


in  the  interval  [0,5ir]  for  the  initial  conditions 
y(0)  =  0 
y’(0)  =  1 

Reduce  the  second  order  linear  differential  equation  to  a  system  of  two  first 
order  equations 

dyl 

IT"  y2 

dy2 

“*T"  -yl 

with  the  initial  conditions 
y1(0)  =  0 
y2(0)  =  1 

Use  the  Adams-Moulton  predictor-corrector  method  to  integrate  the  given 
differential  equations  over  [0,5^].  An  adjustable  step  size  is  used  with  size 
determined  by 

(I)  Relative  truncation  error 

(II)  Absolute  truncation  error 
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O  C 

The  error  window  is  [10  ,10  ]. 

Two  types  of  printing  options  are  illustrated  in  this  example 

(I)  Printing  every  kth  step  (10  in  this  example) 

(II)  Printing  when  the  independent  variable  has  changed  by  a  pre-selected 
amount  (1  in  this  example) 

The  analytic  solution  of  the  given  equation  with  respect  to  the  given 
initial  conditions  is 
y(x)  -  sin(x) 

The  actual  value  of  sin(x)  is  computed  for  each  x  and  compared  to  the  calculated 
value  of  the  solution.  At  every  print  cycle  the  solution,  sin(x),  and  the 
difference  of  the  actual  and  calculated  solutions  are  printed.  These  results 
are  in  Table  2. 

The  main  (calling)  program  and  the  necessary  subroutines  to  accomplish  the 
solution  of  the  stated  problem  in  the  manner  indicated  are  in  Appendix  B. 

Example  No.  3 

Solve  the  system  of  first  order  linear  differential  equations 
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in  the  interval  [0,5it]  for  the  initial  conditions 

y^O)  =  0 

y2(0)  =  1 

y3(0)  -  1 

y4<0)  -  l 

y5(0)  -  0 

Use  the  Adams-Moulton  predictor-corrector  method  to  integrate  the  given 

system  of  differential  equations  over  [0,5tt].  An  adjustable  step  size  is  used 

—8  —5 

with  size  adjusted  to  hold  the  relative  truncation  error  in  the  window  [10  ,10  ]. 

Printing  occurs  every  40  accepted  steps.  The  analytic  solutions  of  the 
given  equations  with  respect  to  the  given  initial  conditions  are 

y-^(x)  =  sin(x) 

y2(x)  =  cos (x) 

y3(x)  =  (x+l)e_x 

y^(x)  =  x+eX 

y5(x)  -  loge(x+l) 

The  actual  values  of  the  solutions  are  computed  at  each  x  and  compared  to 
the  calculated  values  of  the  solutions.  At  each  print  cycle  the  actual  solutions, 
calculated  solutions,  and  differences  are  printed.  These  results  are  in  Table  3. 

The  main  (calling)  program  and  the  necessary  subroutines  to  accomplish 
the  solution  of  the  stated  problem  in  the  manner  indicated  are  in  Appendix  C. 
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VII.  REMARKS 

Some  suggestions  and  warnings  are  in  order  based  on  the  years  of  experience 
with  previous  versions  of  this  routine.  This  version  of  FN0L3  is  a  conversion 
from  the  IBM  7090,  and  what  was  one  auxiliary  routine  with  multiple  entries 
has  been  changed  to  three  auxiliary  routines.  In  addition,  double  precision 
has  been  deemed  unnecessary  due  to  the  increased  precision  of  the  CDC  6400. 

.J:  Adams-Moulton  is  the  prime  method  of  this  routine,  but  Runge-Kutta  is 

needed  to  start,  and  whenever  the  time  step  changes.  Adams-Moulton 
requires  only  2  entries' into  DERIV  per  time  step  while  Runge-Kutta 
requires  4. 

L^:  This  parameter  is  a  hangover  from  previous  versions  of  FN0L3  in  which  the 

user  did  not  write  an  OUTPUT  routine  but  simply  specified  how  many  (L) 
auxiliary  values  were  to  be  printed. 

M:  FN0L3  always  prints  initial  and  final  values  even  if  M  has  been  set  to 

some  very  large  value.  However,  some  output  is  always  needed  to  shed 
light  on  how  the  integration  has  proceeded. 

Printing  at  specified  intervals  cannot  be  insured  except  by  fixed 
timesteps  or  terminating  and  then  printing.  Setting  M  =  0  and 
Y(N+2)=DELTAY  in  MAIN  and  Y(N+1)=Y(1)  in  TERM,  only  causes  no  printing 
until  Y(l)  changes  by  at  least  Y(N+2).  To  force  printing  at  exact 
intervals  of  Y(N+2)  in  Y(l),  put  FN0L3  in  a  loop  in  MAIN  and  in  TERM  set 
T=Y(1)-C.  In  MAIN  set  C=Y(l)+Y(N+2)  before  calling  FN0L3.  C  may  be 
transmitted  via  the  Y  array  or  through  COMMON.  The  final  value  of  Y(l)  of 
each  integration  becomes  the  initial  value  for  the  next  integration. 

Y_:  Unless  COMMON  is  used,  additional  parameters  necessary  for  computation 

should  be  sent  from  routine  to  routine  through  the  Y  array  in  locations 
starting  with  Y(N+4).  If  M  4  0  then  Y(N+1)  and  Y(N+2)  are  available 
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to  the  programmer.  If  XNE  =  0  then  Y(N+3)  is  available.  In  other 
cases  these  three  locations  have  special  meaning.  See  the  discussions 
on  M  and  XNE  above  and  in  Section  V. 

D:  We  repeat  here  for  emphasis.  DO  NOT  USE  a  D  beyond  D(N)!  In  all  we 

have  said,  there  is  an  assumption  that  the  derivatives  can  be  computed 

explicitly  from  X  and  the  current  Y  values.  If  this  is  not  so,  the 

alternative  is  to  perform  root-finding  within  DERIV.  If  the  routine 

that  calls  FN0L3  initializes  the  D  array,  these  values  will  be  available 

in  DERIV.  However,  unless  they  are  saved  in  another  array,  each  time 

DERIV  is  entered,  the  previous  D's  are  lost. 

ERROR:  At  the  start  of  an  integration  or  whenever  the  step  size  changes,  the 

absolute  truncation  errors  are  not  computed  directly  for  the  first 

step.  Runge-Kutta  is  called  at  that  time,  and  its  error  term  requires 

having  yf^,.  and  y^^  ,  see  equation  IV-3.  The  first  term,  y .  , 
i  9 n*  X  i  ^ in i  1  i  ^ rtf- j, 

means  reaching  by  taking  two  steps  of  size  h,  while  the  second 

/  ou  \ 

term,  y^  means  reaching  by  taking  one  step  of  size  2h.  It 

is  only  after  these  three  steps  are  taken  that  an  error  term  can  be 
computed.  Accepting  or  rejecting  the  first  step  is  based  on  this  error, 
which  is  then  treated  as  if  it  were  the  first  error.  At  least  three 
Runge-Kutta  steps  of  size  h  are  taken  before  Adams-Moulton  takes  over. 


VIII.  ACKNOWLEDGMENT 

NOL  has  had  an  ordinary  differential  equations  solver  since  1959,  and 
FN0L3  is  the  current  version.  We  offer  our  sincere  thanks  to  L.  Hieb , 

M.  Vander  Vorst  and  J.  Lyles  who  helped  to  checkout  and  criticise  the  FN0L3 
program  and  writeup.  A  number  of  improvements  were  made  on  the  basis  of  their 
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TABLE  1 


=  uA  cos(a)x) 


Method  (I) 


X 

YT-TRUE 

sin(x) 

YC=C0MPUTED 

sin(x) 

(YT-YC) 

0.0 

0.0 

0.0 

0.0 

0.5 

.479426 

.479426 

-17xl0-9 

1.0 

.841471 

.841471 

-29xl0~9 

1.5 

.997495 

.997495 

-35xl0-9 

2.0 

.909297 

.909297 

-32xl0~9 

2.5 

.598472 

.598472 

-21xl0-9 

3.0 

.141120 

.141120 

-49xl0-10 

3.5 

-.350783 

-.350783 

12xl0-9 

4.0 

-.756802 

-.756803 

26x10  9 

4.5 

-.977530 

-.977530 

34xl0~9 

5.0 

-.958924 

-.958924 

33xl0-9 

5.5 

-.705540 

-.705540 

25xl0-9 

6.0 

-.279415 

-.279416 

97xl0-10 

6.5 

.215120 

.215120 

-75xl0-10 

7.0 

.656987 

.656987 

-23xl0"9 

7.5 

.938000 

.938000 

-33xl0-9 

8.0 

.989358 

.989358 

-34xl0~9 

8.5 

.798487 

.798487 

-28xl0-9 

9.0 

.412118 

.412118 

-14xl0~9 

9.424778 

-39xl0~9 

-4lxl0-9 

17xl0-10 
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TABLE  1 


dx 

cos(o)x) 

Method 

(II) 

YT-TRUE 

YOCOMPUTED 

(YT-YC) 

X 

sin(x) 

sin(x) 

0.0 

0.0 

0.0 

0.0 

0.5 

.479426 

.479426 

22xl0-9 

1.0 

.841471 

.841471 

92xl0~9 

1.5 

.997495 

.997495 

18xl0-8 

2.216974 

.798391 

.798391 

24xl0-8 

2.822212 

.313979 

.313979 

16x10  8 

3.219365 

-.0776944 

-.0776946 

17xl0"8 

3.729709 

-.554795 

-.554795 

13xl0-8 

4.240053 

-.890508 

-.890508 

47xl0-9 

4.750397 

-.999278 

-.999278 

-59xl0-9 

5.260741 

-.853385 

-.853384 

-16xl0"8 

5.771085 

-.490009 

-.490009 

-23xl0~8 

6.220706 

-.0624391 

-.0624388 

-25xl0_8 

6.658950 

.366984 

• 

.366985 

-24xl0-8 

7.155075 

.765546 

.765546 

-18xl0"8 

7.651199 

.979510 

.979510 

-92xlO-9 

8.353455 

.877835 

.877835 

29xl0-8 

9.288410 

. 135945 

.135946 

-6lxl0-8 

9.424778 

-39xl0-9 

62xl0~8 

—66x10  8 
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TABLE  1 


=  uA  cos(wx) 


Method  (III) 


X 

YT=TRUE 

sln(x) 

YCCOMPUTED 

sin(x) 

(YT-YC) 

0.0 

0.0 

0.0 

0.0 

0.5 

.479426 

.479426 

12xl0~9 

1.303726 

.964548 

.964548 

91xl0-9 

2.309936 

.739048 

.739048 

77xl0-9 

3.316146 

-.173668 

-.173668 

-46xl0~8 

4.322356 

-.924896 

-.924895 

-10xl0~7 

5.297378 

-.833718 

-.833717 

-llxlO-7 

6.225618 

-.0575354 

-.0575345 

-83xl0-8 

7.153858 

.764763 

.764763 

-4lxl0-8 

8.082099 

.974094 

.974094 

-23xl0-8 

9.102580 

.316653 

.316653 

-51xl0-8 

9.424778 

-39xl0-9 

74xl0-8 

-78xl0'8 
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TABLE  2 


dx 


2 


-y 


Relative  Truncation  Error 


X 

YC=C0MPUTED 

sin(x) 

YT'TRUE 

sin(x) 

(YT-YC) 

0.0 

0.0 

0.0 

0.0 

0.381345 

.372169 

.372169 

27xl0~9 

1.194792 

.930140 

.930139 

-22xl0~8 

2.008240 

.905838 

.905838 

-52xl0-8 

2.821688 

.314477 

.314476 

-39xl0"8 

3.322058 

-.179488 

-.179488 

llxlO-9 

4.055814 

-.792088 

-.792087 

72xl0-8 

4.710873 

-1.00000 

-.999999 

12xlO-7 

5.262094 

-.852680 

-.852679 

llxlO-7 

5.953946 

-.323324 

-.323324 

61x10  8 

6.645797 

.354718 

.354717 

-30xl0-8 

7.337648 

.869636 

.869635 

-12x1 0-7 

8.029500 

.984638 

.984636 

-17xlO~7 

8.721351 

.646836 

.646835 

-13xlO-7 

9.413203 

.0115755 

.0115752 

-32xl0-8 

10.105054 

-.629009 

-.629008 

99xl0-8 

10.796905 

-.980332 

-.980330 

20xl0~7 

11.488757 

-.880833 

-.880831 

21xl0-7 

12.180608 

-.376267 

-.376266 

12xl0~7 

12.872459 

.301332 

.301332 

-34xl0-8 

13.564311 

.840358 

.840356 

-19xl0-7 

14.256162 

.992931 

.992928 

-27xl0-7 

14.948014 

.688887 

.688885 

-22xlO“7 

15.639865 

.0680464 

.0680457 

-69xl0-8 

15.707963 

76xl0-8 

27xl0-8 

-50xl0-8 

21 


NOLTR  71-2 


TABLE  2 


Absolute 

Truncation  Error 

X 

YC-COMPUTED 

sin(x) 

YT=TRUE 

sin(x) 

(YT-YC) 

0.0 

0.0 

0.0 

0.0 

1.008358 

.845958 

.845957 

-91xl0~8 

2.030260 

.896293 

.896290 

-30xl0-7 

3.052163 

.0893123 

.0893108 

-16xl0~7 

4.074065 

-.803100 

-.803096 

46xl0~7 

5.095968 

-.927340 

-.927331 

89xl0-7 

6.117870 

-.164568 

-.164563 

43xl0-7 

7.139773 

.755619 

.755612 

-75xlO-7 

8.161675 

.953050 

.953035 

-15xl0-6 

9.183578 

.238876 

.238868 

-79xl0-7 

10.205480 

-.703788 

-.703778 

96xl0-7 

11.227383 

-.973273 

-.973253 

21xl0-6 

12.249285 

-.311811 

-.311799 

12x10  8 

13.271187 

.647905 

.647894 

-llxlO"6 

14.293090 

.987895 

.987869 

-27xl0-6 

15.314992 

.382952 

.382935 

-17xl0-6 

15.707963 

76xlO-7 

27xl0-8 

-74xl0-7 
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TABLE  3 


X 

EQUATION  NO. 

YC= COMPUTED  y 

YT=TRUE  y 

(YT-YC) 

0.0 

1 

0.0 

0.0 

0.0 

2 

1.0 

1.0 

0.0 

3 

1.0 

1.0 

0.0 

A 

1.0 

1.0 

0.0 

5 

0.0 

0.0 

0.0 

2 . 3055 

1 

. 7A2059 

. 7A2058 

-13xl0-8 

2 

-.670335 

-.670335 

91xl0-9 

3 

.329598 

.329598 

-20xl0~9 

A 

12.33A2 

12.33A2 

15xl0~7 

5 

1.19557 

1.19557 

19xlO-8 

A.  6901 

1 

-.999753 

-.999753 

3AxlO~8 

2 

-.0222392 

-.0222391 

63xl0-9 

3 

-.0522659 

-.0522659 

-37xlO~10 

A 

113.559 

113.559 

35xl0'6 

5 

1.7387A 

1.7387A 

20xl0-8 

7.07A8 

1 

.711517 

.711516 

-31xl0-8 

2 

.702670 

.702669 

-A3xl0-8 

3 

.00683235 

.00683235 

-18xl0-11 

-5 

A 

1188.93 

1188.93 

59x10 

-8 

5 

2.08875 

2.08875 

20x10 
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X 

EQUATION  NO. 

YC=COMPUTED  y 

YT-TRUE  y 

(YT-YC) 

9.4595 

1 

-.0347456 

-.0347457 

-98xl0"9 

2 

-.999397 

-.999396 

7lxl0"8 

3 

.815249xl0~3 

.815249xl0-3 

43xl0-12 

4 

12839.3 

12839.3 

86xl0-4 

5 

2.34751 

2.34751 

20xl0-8 

11.8442 

1 

-.660999 

-.660999 

70xl0~8 

2 

.750388 

.750387 

-56x10  8 

3 

.922205xl0-4 

.922205xl0~4 

-12 

14x10  x 

4 

139289. 

139289. 

12xl0-2 

5 

2.55289 

2.55289 

20xl0-8 

14.2289 

1 

.995795 

.995794 

-llxlO-7 

2 

-.0916176 

-.0916177 

-90xl0"9 

3 

.100723xl0-4 

.100723xl0-4 

27xl0~13 

4 

.151197xl07 

. 151197xl07 

15xl0_1 

5 

2.72320 

2.72320 

20xl0-8 

15.7080 

1 

.7977l8xl0-6 

.587949xl0~6 

-21xlO~8 

2 

-1.00000 

-1.00000 

12xl0-7 

3 

. 251792xl0-5 

. 251792x10” 3 

87xl0~14 

4 

.663563xl07 

.663564xl07 

75xlO-1 

5 

2.81589 

2.81589 

20xl0-8 
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APPENDIX  A 

LISTING  OF  EXAMPLE  NO.  1  WITH  CONTROL  CARDS 
FOR  THE  NOLOS  SYSTEM  USED  AT  NOL 


CCAFNOL, T 100, CM60000. SYSTEM, 039, ORLOW. 

ATTACH! ABC, NOLDIN) 

COPYN(0»DEF,ABC)  FIRST  PARAMETER  IS  THE  NUMBER  »0' 

RETURN! ABC) 

FTN! L ) 

LOAD ! LGO ) 

DFF  . 

'  RECORD  SEPARATOR  c  (7-8-9)  PUNCH  IN  COL.  1 

REWIND! ABC) 

FN0L3 , 1 , ABC 

•  RECORD  SEPARATOR  =  (7-8-9)  PUNCH  IN  COL.  1 

PROGRAM  EX  1  (  INPUT , OUTPUT ) 

C 

C  FXAMPLF  NO.  1 

C  SOLVE  THE  INITIAL  VALUE  PROBLEM  DY/Dx=A*W*COS ( W*X ) 

C  USING  THE  FOLLOWING  METHODS 

C 

C  (I)  4TH  ORDER  RUNGE-KUTTA  (CONSTANT  STEP) 

C  (II)  4TH  ORDtR  ADAM5-M0ULT0N  METHOD  (VARIABLE  STEP) 

C  (III)  4TH  ORDER  RUNGE-KUTTA  (VARIABLE  STEP) 

C 

DIMENSION  Y ( 20 )  , D  <  2  0 ) 

EXTERNAL  DERIV, TERM, OUTPUT 
COMMON  A , W 
A  =  1  . 

W=  ]  . 

DO  200  J  =  1 »  3 

C  NUMBER  OF  FOUATIONS 

N=  1 

C  INITIAL  STFP  SIZE 

G=  •  1 

C  NUMBER  OF  EXTRA  Y'S 

L  =  2 

C  NUMBER  OF  ACCEPTED  STEPS  TAKEN  BETWEEN  PRINTS 

M  =  5 

C  FRROR  WINDOW  EXPONENT 

<NE=5. 

Y (  1  )  =0 . 

Y ( N+3 )  =  • 0  1 
X  —  0  • 

PRINT  1000 

1000  FORMAT  ( 1 HI , 4X , 1 HX , 1 7X , 9H YT  =  7  RUE  Y , 1 2  X ♦ 1 3HYC  =  COMPUT  ED  Y»11X, 
1  10H  (YT-YC) ) 

CALL  FN0L3  ( J , N ,G , L ,M , XNE , X , Y * D , DER I  V , T ERM , OUTPUT ) 

200  CONTINUE 
STOP 
FND 
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SUBROUTINE  DERIV  (X.Y.D) 

DIMENSION  D ( 1 ) 

COMMON  A  » W 
D (  l)»A*W*COS(W*X) 

RETURN 

END 

SUBROUTINE  TERM  (X.Y.D.T) 

DTMFNSION  Y(l) 

COMMON  A  » W 
T=X-9. 424778 

Y  < 5 ) =A*S I N ( W*X ) 

Y  (  6  )  *Y  (  5  )  — Y  (  1  ) 

RETURN 

END 

SUBROUTINE  OUTPUT  ( X , Y . D , ERROR , N , L , H ) 

DIMFNSION  Y ( 1 ) 

PRINT  1000.  X»Y(5)*Y(1)»Y(6) 

1000  FORMAT  ( 1H  » FI  0 .6  •  1  OX . 3 ( E13 • 6 » 10X ) ) 

RETURN 

END 

ORLOW  039  CCAFNOL  END  OF  FILE  =  16-7-8-9)  PUNCH  IN  COL. 
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LISTING  OF  EXAMPLE  NO.  2 


PROGRAM  EX2 ( INPUT  » OUT PUT ) 

C  FXAMPLE  NO.  2 

c 

C  SOLVE  D**2Y/DX**2=-Y 

C 

C  USING  ADAMS-MOULTON  PREDICTOR  CORRECTOR  WITH 

C  STEP  SIZE  CONTROLLED  BY 

C 

C  (  I  )  RELATIVE  ERROR 

C  (II)  ABSOLUTE  ERROR 

C 

DIMENSION  Y ( 30 ) .D(30> 

EXTERNAL  DER I  V , TERM .OUTPUT 
C  ADAM S-MOUL TON  METHOD 

J  =  2 

C  NUMBER  OF  EQUATIONS 

N*  2 

C  MINIMUM  RELATIVE  ERROR  DIVISOR 

Y ( N+3 ) = . 00 1 

C  NUMBER  OF  EXTRA  Y'S 

L  =  2 

C  NUMBER  OF  ACCEPTED  STEPS  bETWELN  PRINT  CYCLES 

M*  1 0 

C  EXPONENT  ERROR  WINDOW 

XNE  =  5 . 

C  INITIAL  STEP  SIZE 

G=  .01 

C  INITIAL  CONDITIONS 

X  =  0. 

Y (  1  )  =0. 

Y ( 2 ) =1  . 

PRINT  1000 

1000  FORMAT  (1H1»4X»1HX*13X»1 3HYC “COMPUTED  Y»16X.RHYT=TRUE  Y»11X» 
1  10H  (YT-YC)) 

CALL  FN0L3(J»N»G»L»M»XNE*X*Y»D»DERIV»TERM .OUTPUT ) 

C  RESET  THE  PRINTING  OPTION 

M  =  0 

Y ( N+2 )  =  1  • 

C  RESET  THE  INITIAL  CONDITIONS 

X  =  0. 

Y ( 1 ) =0 . 

Y  (  2  )  *  1  • 

C  FHANGF  TO  ARSOLUTF  ERROR 

Y ( N  +  3  >  =-Y ( N  +  3 ) 

PRINT  1000 

CALL  FN0L3 (J»N»G»L»M*XNE»X»Y,D»DERIV. TERM .OUTPUT ) 

STOP 

END 
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SUBROUTINE  DERIVIX.Y.D) 

DIMENSION  Y  < 1 ) » D ( 1 ) 

D( 1 ) =Y( 2  ) 

D ( 2 ) =-Y (  1  ) 

RETURN 

END 

SUBROUTINE  TERM ( X , Y » D , T ) 

DIMENSION  Y ( 1 ) 

T  =  1 5 . 70  7963-X 

Y  (  3  )  =X 

Y  C  6 ) =S I N ( X ) 

Y  (  7 ) =Y  <  6 ) -Y (  1  ) 

RETURN 

END 

SUBROUTINE  OUT PUT ( X , Y  » D , ERROR ,N,L,H) 
DIMENSION  Y ( 1 ) 

LF=N+4 
LL  =N+3  +  L 

PRINT  1000»  X, Y( 1 ) , ( Y( I ) ♦ I=LF ,LL) 
1000  FORMAT ( 1H  *  FI  0*6 • 10X t 3 ( El  3 .6 » 10X ) ) 
RETURN 
END 
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PROGRAM  EX3 ( INPUT .OUTPUT ) 

C  EXAMPLE  NO.  3 

C 

C  SOLVE  THE  SYSTEM  OF  FIRST  ORDER  EQUATIONS 

C 

C  DY ( 1 ) /DX  =  Y ( 2  ) 

C  DY ( 2 ) /DX=-Y ( 1 ) 

C  DY ( 3 ) /DX* ( 1 ./( X+l. )-l. >*Y< 3 ) 

C  DYI4)/DX=Y(4)-X+1. 

C  DY (5) /DX=1./(X+1. ) 

C 

C  USING  ADAMS-MOULTON  PREDICTOR  CORRECTOR  METHOD  WITH 

C  STEP  SIZE  CONTROLLED  BY  RELATIVE  ERROR 

C 

DIMENSION  Y ( 30 ) ,D(30) 

EXTERNAL  DER IV , TERM .OUTPUT 
C  ADAMS-MOULTON  METHOD 

J~2 

C  NUMBER  OF  EQUATIONS 

N  =  5 

C  MINIMUM  RELATIVE  ERROR  DIVISOR 

Y  (  N  +  3 )  =  . 0 1 

C  NUMBER  OF  EXTRA  Y'S 

L=  1 0 

C  NUMBER  OF  ACCEPTED  STEPS  BETWEEN  PRINT  CYCLES 

M=40 

C  EXPONENT  ERROR  WINDOW 

XNE  =  5 • 

C  INITIAL  STEP  SIZE 

G=  •  02 

C  INITIAL  CONDITIONS 

X  =  0. 

Y ( 1 ) =0 . 

Y  (  2  )  =  1 . 

Y  (  3  )  =  1  . 

Y  (  4  )  =  1  . 

Y( 5)=P. 

PRINT  lOf'O 

100^  FORMAT (  1H1  *4X*1HI  *13X*14HYC=C OMPUT FD  Y  I  ,  1 4X , 1 OHYT  =  TRUE  YI,11X, 
1  10H  (YT-YC)) 

CALL  FN0L3IJ  *N*G,L »M, XNE, X»Y,D,DER IV, TERM, OUTPUT) 

STOP 

END 
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SUBROUTINE  DERIV(X,Y,D) 

DIMENSION  Y ( 1 ) * D ( 1 > 

Dd  )  =  Y(  2  > 

D ( 2 ) =-Y ( 1 ) 

D(3)  =  d./(X+1.)-1.)*Y<3) 
D(4)*Y(4)-X+1. 

D ( 5 ) = 1 • / ( X+l  .  ) 

RETURN 

END 

SUBROUTINE  T ERM ( X , Y , D , T ) 

DIMENSION  Y ( 1 ) 

T=X-5. *3.141892536 
Y ( 9 ) =S I N ( X ) 

Y ( 10I=COS(X) 

Y (  11)  =  (X+1. )  *EXP ( -X ) 

Y ( 12 ) =X+EXP ( X ) 

Y ( 1 3 ) =  A  LOG ( X  +  l . ) 

DO  10  1=14,18 

Y(  I  )=Y(  I -5  >  — Y ( 1-13) 

10  CONTINUE 
1ETURN 
END 

SUBROUTINE  OUTPUT ( X ,Y ,D .ERROR ,N ,L ,H) 
DIMENSION  Y(l) 

PRINT  1000,  X 

1000  FORMAT ( 1H0 ,50X ,3HX=  ,F7.4) 

J=  N  +  3 
tc  =  ?#N+? 

DO  20  1=1, N 
J  =  J+1 
K  =  K+ 1 

PRINT  2000,  I  ,  Y (  I  ) , Y ( J ) , Y ( K ) 

2000  FORMAT ( 1H  ,  I  5 » 1 5X  ,  3 ( E 1 3 . 6 , 1 OX ) ) 

20  CONTINUE 
RETURN 
END 
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SUBROUTINE  FNOL3 ( J , NN , G ,L ,MPR , XNC » X , Y ,D ,DER I  V . TERM, OUTPUT ) 

C  001  ^INTEGRATION  METHOD  CONTROL 

C  002  NN-NUMBER  OF  SIMULTANEOUS  DIFFERENTIAL  EQUATIONS 
C  003  G=F I RST  INTERVAL  OF  INTEGRATION 
C  004  L-NUMBER  OF  YSS  GREATER  THAN  N  TO  BE  PRINTED 
C  005  MPR=PR I  NT  FREQUENCY 

C  006  XNE-CONTROL  FOR  INTERVAL  OF  INTEGRATION 
C  007  X* I NDEPENDENT  VARIABLE 
C  008  Y*DEPENDENT  VARIABLE 
C  009  D« ARRAY  CONTAINING  DERIVATIVES 

C  010  DERI  V* SUBROUT  I NE  IN  WHICH  DIFFERENTIAL  EQUATIONS  ARE  FOUND 
C  Oil  TERM<=  SUB  ROUT  I  NE  FOR  TERMINATION  CONDITION 
C  012  OU  TPUT  =  SUB ROUT  I NE  FOR  PRINTING  OUTPUT 
C  NRET=TERM  LOOP  COUNTER 

C  NPR=COUNT  OF  STEPS  SINCE  LAST  PRINT 

C  PC* Y ( N+ 1 )= PRINT  CONTROL  OTHER  THAN  STEP  COUNT 

C  JJ* J-2  =  0  FOR  AM ♦ - 1  FOR  RK,  +  1  FOR  RK2 

C  MK*AM  RK  STEP  COUNT 

C  XD  =  DP  FORM  OF  X 

C  THIS  SUBROUTINE  CAN  ALSO  BE  RUN  IN  DOUBLE  PRECISION  BY  REMOVING 
C  THE  *C»  IN  COLUMN  1  ON  THE  DOUBLE  PRECISION  STATEMENT 
C  DOUBLE  PRECISION  XD , YD , YP , YC , YK , H , HC , XS . XD5 . HS »HN ,HD24 ,HD6 

D I  MENS I  ON  C( 3) »Y( 30 ) » YD( 30)  »YP  C  30) t YC ( 30  >  *D( 50) *DM( 30»4) »DK(30*4) 
1 »  ERROR ( 30) »YK( 30) 

DATA  EP6,EP11»M4/1.E-6,1.E-11 .-4/ 

DATA  <C(K)  ,K=1 ,3) /2*.5,1 ./ 

DATA  HMAX5/1 .F35/ 

NHT  S  =  0 
FP2=0. 

N*NN 
J J= J-2 
H=G 
HN  =  H 
MK  =  1 
NRET*M4 
JTEST= 1 

IF  (JJ  .LT.  0)  JTEST=4 
IF  (XNE.EQ.  0.)  GO  TO  15 
EC=Y ( N+3 ) 

EUP=10.**(-XNE) 

EL0=EUP*.00l 
EM«EL0»31. 6227766 
15  XD=X 
XS  =  XD 

DO  2  0  1=1  *  N 
ERROR ( I ) =0 • 

20  YD (  I »  = Y ( I ) 

CALL  DERIV(X*Y»D) 

CALL  TERM  (X»Y»D.F) 

C  PRINT 
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50  CALL  OUTPUT! X»Y»D  » ERROR  , N  ,  L  ,H  ) 

IF  (NRET)  65  »60 ,55 
55  PRINT  3000*  HN 

3000  FORMAT ( 108H1 EXECUT ION  TERMINATED  BECAUSE  INTERVAL  OF  INTEGRATION  L 
1  ESS  THAN  l.OE  -6  TIMES  INDEPENDENT  VARIABLE  (X).  H  =,1PE15.7) 

STOP 

60  RETURN 
65  NPR  =  0 

IF  ( MPR  .LE.  0)  PC=  Y ( N+l ) 

100  IF  ( JTCST  .EO.  5  .AND.  H  .NE.  HN)  GO  TO  455 
IF  (JJ  .GE.  0)  H=HN 

IF  (MK  .NE.  0  • OR .  JJ  .NE.  0)  GO  TO  300 

C - THE  ADAMS  MOULTON  METHOD 

200  HD24=H/24. 

JAM  =  0 

DO  210  1=1 ,N 

YP I = ( 55 • *DM ( I , 1 ) +37 • *DM ( 1,2) ) -( 59.*DM< 1,3 )+9.*DM( 1,4)) 

YP ( I ) =YD ( I ) +HD24* YP I 

Y  (  I  )  = YP ( I ) 

210  CONTINUE 

X=XD+H 

CALL  DERI  V! X ,Y ,DM( 1 ,4 )  ) 

DO  220  1  =  1  , N 

YP I = ( 9 . *DM ( I ,4)+19.*DM( I , 1 )+DM< 1,2) )-5.*DM( 1,3) 

YC( I ) = YD ( I ) +HD24*YP I 
ERROR! I ) = ( YP ( I ) -YC ( I ) )/14. 

C  THIS  ADDS  IN  A  2D  CORRECTION 
YC! I ) =YC( I ) +ERROR < I ) 

220  CONTINUE 

IF  (XNE.NE..0)  GO  TO  700 
GO  TO  455 

C - THE  RUNGE  KUTTA  METHOD 

300  GO  TO  (301,309,308,309,303)  ,JTEST 

301  DO  302  1  =  1  , N 
YK( I ) = YD ( I ) 

302  CONTINUE 
XDS=  XD 

GO  TO  309 

303  DO  304  1  =  1  , N 
YK( I ) = YC ( I  ) 

304  CONTINUE 
XDS=XD+H 

308  HS=H 
H=2.*H 

GO  TO  320 

309  X=XD 
J  A  M  =  1 

DO  310  1=1, N 

Y  (  I  )  =YD ( I ) 

DK( I ,1 ) =D( I  ) 

310  CONTINUE 
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IF  (JTEST  .LE.  2)  CALL  DER I  V ( X , Y , DK ) 

IF  (MK  .GT.  1  .OR.  JTEST  .GT.  1)  GO  TO  320 
DO  315  1  =  1  ,N 
DM ( I ,4)=DK( I  .1  ) 

315  CONTINUE 
320  JO  335  K  =  2  *  4 
HC=H*C( K-l ) 

DO  330  I  =  1  ♦  N 

Y (  I ) =YD ( I )  +  HC*DK ( I  *  K-l ) 

330  CONTINUE 
X  =  X  D+HC 

CALL  DE RIVfX.Y* DK ( 1 » K )  ) 

335  CONTINUE 
HD6*H/6 • 

DO  340  1  =  1  , N 

YPI=DK(  I  ,1 ) +  DK ( I » 4 ) +2 • * ( DK ( I , 2 )+DK( 1,3)) 
YC( I ) *YD( I ) +HD6*YPI 
340  CONTINUE 

GO  TO  ( 360,390,370,455,370) , JTEST 
360  DO  365  1  =  1  , N 
YP ( I ) = YC ( I ) 

365  CONTINUE 
JTEST*3 
GO  TO  308 
370  DO  380  1  =  1  , N 
YD ( I ) =YP ( I ) 

YP (  I  ) =YC( I  ) 

380  CONTINUE 
H*HS 
XD=XD+H 
JT  EST  =  2 

IF  (MK  .EQ.  1)  GO  TO  309 
GO  TO  451 
390  DO  400  1  =  1  » N 

ERROR ( I ) = { YC ( I ) — YP ( I ) ) / 1 5 . 

YC(  I )=YC( I  ) +ERROR {  I  ) 

YP ( I ) =YC ( I  ) 

400  CONTINUE 
JTEST=5 

IF  (XNE.NE..O)  GO  TO  700 
C - ACCEPT  THE  STEP  SIZE 

450  IF  (JAM  .EQ.  0)  GO  TO  455 

IF  (MK  .EQ.  3  .AND.  JJ  .EQ.  0)  GO  TO  455 
IF  (MK  .NE.  1)  GO  TO  303 
IF  (JTEST  .EQ.  1 )  GO  TO  455 

451  DO  452  1  =  1  , N 
Y (  I ) =YD  <  I  ) 

452  CONTINUE 
GO  TO  466 

455  DO  459  NQ=1  ,N 
YD ( NO ) =  YC ( NQ ) 
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Y ( NQ )  =  YD  ( NO ) 

459  CONTINUE 

IF  (JJ  .GE.  0)  JTf ST= 1 

IF  (MK  .NE.  0  .OR.  JJ  .NE.  0  -OR.  NRET  .NE.  M4 )  GO  TO  465 
DO  460  1  =  1  » N 
DM { I  ,4) =DM< I * 2 ) 

DM (  1*2)*  DM  C 1*3) 

DM ( I »  3 ) =DM ( I  ,1 ) 

460  CONTINUE 

465  XD=XD+H 

466  X=XD 

IF  (MK  .EQ.  3)  MK  =  0 
CALL  DERIV(X,Y,D) 

DO  470  1=1, N 
DM (  I  * MK  +  1 ) =D(  I  ) 

470  CONTINUE 
FP  =  F 

CALL  TERM  (X,Y,D,F) 

C - DO  YOU  TERMINATE 

500  IF  (ABS(F)  .LE.  EP6 )  GO  TO  800 
IF  (FP  .EQ.  0.)  GO  TO  510 

IF  (NRET  .NE.  M4  .OR.  F*FP  .LT.  EP11)  GO  TO  805 
510  XS=XD 

IF  (MK  .NE.  0  .AND.  H  . EO .  HN )  MK=MK+1 

C - DO  YOU  PRINT 

600  NPR  =  N'PR  +  1 

IF  ( MPR  .LE.  0)  GO  TO  610 
IF  ( NPR  .GE.  MPR)  GO  TO  50 
GO  TO  100 

610  IF  ( ABS ( Y( N+l ) -PC )  .GE.  Y(N+2))  GO  TO  50 
GO  TO  100 

C - DETERMINING  THE  STEP  SIZE 

700  HB  =  HMAX5 

DO  760  I  =  1  *N 
Z  =  ABS ( E RROR  (  I  )  ) 

IF  (YC(I)  .EQ.  0.)  GO  TO  720 
ZZ  =  YC (  I  ) 

ZZ=ABS ( ZZ ) 

IF  (EC)  720*710*705 
705  IF  (EC  .GT.  ZZ)  ZZ=EC 
710  Z= Z/ZZ 

720  IF(Z.GT.ELO.AND.Z.LT.EUP)  GOTO  750 
HB  =  AMIN1 (HB,EM/(Z+EPU ) ) 

GO TO 760 

750  HB=AMI N1 (HB, 1 . ) 

760  CONTINUE 

IF  (HB  .NE.  1. )  GO  TO  765 
NHT  S=0 
GO  TO  450 
765  HN=H*H0**. 2 

IF  (MK  .NE.  1)  JT  EST  =  1 
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MK  =  1 

IF  (HB.LT. 1 . )  GOTO  770 
IF  (ABS(HN)  .GT.  ABS(4.*H))  HN=4.*H 
NH  TS=0 
GOTO  450 

770  HEPS=ABS( X*EP6 )  +  EP1 1 

IF  (ABS(HN)  .LT.  A  B  S ( H / 4 . ) )  HN=H/4. 

IF  (ABS(HN)  .GT.  HEPS)  GO  TO  790 

NHTS  =  NHTS  +  1 

IF  (NHTS  .LE.  10)  GO  TO  780 

NRET  =  1 

GO  TO  450 

780  HN=SIGN(HEPS»HN) 

IF  (NHTS  .GT.  1 )  GO  TO  450 
*^oo  IF  (NHTS  .GT.  I)  NHTS*0 

IF  (JAM  .FO.  0)  GO  TO  100 
DO  795  1  =  1  , N 
YD(  I  )  = YK ( I  ) 

*'95  CONTINUE 
XD=XDS 
JTEST= 1 
GO  TO  100 

C - THE  TERMINATION  LOOP 

800  NRET=0 

805  IF  (NRET  .LT.  0)  GO  TO  806 
H=XD— XS 

GO  TO  50 

806  IF(F*FP.LT.0. )  GOTO  810 
IF  (F*FP2.LT.O. )  GOTO  820 
GO  TO  800 

810  FP2  =FP 
HP  =H 
GOTO  830 
820  FP  =FP2 

HP  =H  +  HP 
830  NRE  T  =  NRET+  1 

H=HP*F/ ( FP-F ) 

JTEST  =  4 
GOTO  300 
END 
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FNOL3  VS.  FNOL2 


2-1.  FNOL3  does  not  allow  the  step  size  to  be  changed  by  more  than  a  factor 

of  four.  In  FN0L2,  the  step  size  is  allowed  to  change  by  any  factor,  which 
can  cause  the  step  size  to  be  too  large  or  too  small  in  some  cases. 

E-2.  If  the  step  size  is  too  small,  FN0L3  accepts  the  step  and  the  step  size 
is  increased  for  the  next  step.  FN0L2  does  not  accept  the  step. 

E-3.  FN0L3  calls  the  output  routine  only  if  the  step  was  accepted.  FN0L2 
calls  the  output  routine  and  then  checks  to  see  if  the  step  is  acceptable. 

E-4.  In  FN0L3,  the  truncation  errors  are  added  at  each  step  to  the  dependent 
variables.  FN0L2  does  not  add  the  truncation  errors  to  the  dependent  variables. 

K-5.  FN0L3  allows  the  step  size  to  be  changed  when  using  the  Runge-Kutta 
method.  The  step  size  cannot  be  changed  using  the  Runge-Kutta  method  in  FN0L2. 

K-6.  In  FN0L3,  the  option  was  added  that  if  relative  error  is  used  to  adjust 
the  step  size  and  the  user  sets  Y(N+3)*e(e  >  0),  then  whenever  | Y (I) |  <  Y(N+3), 
relative  error  is  calculated  by  dividing  truncation  error  by  Y(N+3)  rather  than 
Y(I) .  Otherwise  perfect  accuracy  would  be  required  of  every  variable  passing 
through  0. 

F.-7.  FN0L3  uses  XNE,  a  floating  point  rather  than  an  integer,  to  allow  more 
control  over  the  step  size. 
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FNOL3  ARGUMENT  LIST  SUMMARY 


ARGUMENT 

TYPE 

INPUT 

OUTPUT 

REMARKS 

J 

INTEGER 

YES 

The  integration  method 

N 

INTEGER 

YES 

The  number  of  equations 

G 

REAL 

YES 

The  initial  step  size 

L 

INTEGER 

YES 

The  additional  output  Y's 

M 

INTEGER 

YES 

The  print  frequency 

XNE 

REAL 

YES 

The  control  for  step  size 

X 

REAL 

YES 

YES 

The  independent  variable 

Y 

REAL 

YES 

YES 

The  dependent  variables 

D 

REAL 

NO 

YES 

The  derivatives 

DERIV 

TERM 

OUTPUT 

— 

CALL  DERIV (X,Y,D) 

CALL  TERM(X,Y,D,T) 

CALL  OUTPUT (X,Y,D, ERROR, N,L,H) 

F-l 
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